# -*- cording:UTF-8  -*-

# from Standard library
from math import log10,ceil, floor
from sys import exit
import threading, math

# from Numpy & Scipy & pylab
from numpy import *
from scipy import *
from scipy import linalg,special,integrate  
from pylab import find, figure, show

class node(object):
   def __init__(self,nodeX,nodeY):
       self.x = nodeX
       self.y = nodeY

class line(object):
    def __init__(self,node1x,node1y,node2x,node2y, admit):
        self.node1=array([float(node1x),float(node1y)])
        self.node2=array([float(node2x),float(node2y)])
        self.admit=admit
        self.vec=self.node2-self.node1

# Boundary Element & Boundary Element Parameter
class element(object):
    # constructor
    def __init__(self,node1x,node1y,node2x,node2y,admit):
        # place information
        self.node1=array([node1x,node1y])
        self.node2=array([node2x,node2y])
        self.admit=admit

    # calcration each Parameter of each element parameter
    def calcPara(self):
        # objective information
        self.mid=(self.node2+self.node1)*0.5 # middlePoint
        self.range=(self.node2-self.node1)*0.5
        self.length=linalg.norm(self.node2-self.node1)
        self.normalVec=array([self.range[1], -self.range[0]])*2.0/self.length
    
    def __str__(self):
        strings = """element info:
        node1:  %s
        node2:  %s
        admit: %s    
        """ % (self.node1, self.node2, self.admit)
        return strings

class Scatterer(object):
    """
    """
    def __init__(self, filename, admit=0, edge=[]):
        self.set(filename, admit, edge)
    
    def set(self, filename, admit, edge):
        self.filename = filename
        self.lines = self.makeLines(admit, edge)
    
    def setAdmit(self, targetLine, admit):        
        line = self.findLine(targetLine)
        line.admit = admit
    
    def makeLines(self, admit=0, edge=[]):
        nodeList = self.makeNodes(self.filename)
        l = len(nodeList)
        if edge == []:
            temp1 = range(0, l)
            temp2 = range(1, l)
            temp2.append(0)
            edge = array([temp1, temp2]).T
        lines = []
        for ind in edge:
            lines.append(line(nodeList[ind[0]].x, nodeList[ind[0]].y, nodeList[ind[1]].x, nodeList[ind[1]].y, admit))
        return lines
    
    def makeNodes(self, filename):
        lines = self.readNodes(filename)
        nodeList = []
        for line in lines:
            line = line.split(",") 
            nodeList.append(node(float(line[0]), float(line[1])))
        return nodeList
    
    def readNodes(self, filename):
        file_object = open(filename)
        try:
            lines = file_object.read().splitlines()
        finally:
            file_object.close()
        lines = [x for x in lines if x != ""]
        return lines 
    
    def findLine(self, targetLine):
        for line in self.lines:
            if ((line.node1==targetLine.node1)*(line.node2==targetLine.node2)).all():
                return line
        raise ValueError, "Can't find this line."
    
    def plotLines(self, srcPos=[], recPos=[]):
        fig = figure()
        ax = fig.add_subplot(111, autoscale_on=False)
        #ax.hold(True)
        xLim = [0,0]
        yLim = [0,0]   
        for line in self.lines:
            minX = min(line.node1[0], line.node2[0])
            minY = min(line.node1[1], line.node2[1])
            maxX = max(line.node1[0], line.node2[0])
            maxY = max(line.node1[1], line.node2[1])
            
            if xLim[0] > minX:
                xLim[0] =minX
            if xLim[1] < maxX:
                xLim[1] =maxX
            if yLim[0] > minY:
                yLim[0] =minY
            if yLim[1] < maxY:
                yLim[1] =maxY
    
            vec = (line.node2 - line.node1)#*(1/abs(line.node2 - line.node1))
            middle = 0.5*(line.node2+line.node1)
            vec = array(([-vec[1], vec[0]])) # rotation vec 90 degree.
            vec = vec + middle
            
            p = vstack([line.node1, line.node2])
            if line.admit==0:
                ax.plot(p[:,0], p[:,1], "b",label="duct", linewidth=1.5)
                #ax.annotate("", middle, xytext=vec, arrowprops=dict(shrink=1))
            else:
                ax.plot(p[:,0], p[:,1], "r", label="duct", linewidth=2)
                
        if [] not in [srcPos, recPos]:
            ax.plot([srcPos[0]], [srcPos[1]], "or", label="sourcePos")
            ax.plot(recPos[:,0], recPos[:,1], "pb", label="recPos")
        
        plotRange = max(xLim[1]-xLim[0], yLim[1]-yLim[0])
        ax.set_xlim(-plotRange*.1, plotRange*1.1)
        ax.set_ylim(-plotRange*0.5, plotRange*0.5)
        #ax.axis("equal")
        ax.grid(True)
        
            


class BEM(object):
    """
    NAME 
        BEM
    DISCRIPTION
        Conduct BEM calculation for single sound field,  
        where normalized SPL is obtained.
    
    ==============================================================================        
    
    input>>>
    srcPos : position of a line source
    recPos : positions of receivers
    line : lines defining boundaries
    freq : List of frequency Num for calculation.
    temp : temperature of the air.
    limElmNum : maximum number of elements 
    ndiv : (element length) <= (wave length) / ndiv
     
    output>>>
    BEM.result : calculated normalized SPL (freq x recPos)
    
    ==============================================================================        
     
    Based mFile is written by : ISHIZUKA Takashi
    Converted to Python by : ABE Louis 
    Last modified : 2009/01/31 
        
    ==============================================================================        
    
    class-method tree:
    
     BEM - BemThread -+- discBoundElm
                      |
                      +- sbemEqn -+- elmIntGen --- integrand 
                      |           |
                      |           +- elmIntRonElm --- singularIntH02
                      |           |         
                      |           +- greenFunc
                      | 
                      +- calField -+- elmIntGen --- integrand 
                                   |
                                   +- greenFunc

     Ref[1] : Tomonao Okubo, Study on the barrier with acoustically soft    
                 cylinder, PHD thesis ( Kyushu Institute of Design ),      
                 pp. 41-46, 1998.                                          
                  ( on the limit of element length )
    """
    def __init__(self, srcPos, recPos, lines, freqs, temp=15, limElmNum=4000, ndiv=8):
        self.set(srcPos, recPos, lines, freqs, temp, limElmNum, ndiv)
        
    def set(self, srcPos, recPos, lines, freqs, temp, limElmNum, ndiv):
        attributesFromDict(locals())
        self.c = 331.5*sqrt(1+temp/273.0) # Sound Speed
        self.result = zeros([len(self.freqs), self.recPos.shape[0]])
        # check minimum Line length
        self.minLineLen = 1.0*10**(6)  # minimum line lenngth
        for line in self.lines:
            line.length = linalg.norm(line.node2-line.node1);
            if line.length < 1.0*10**(-12): 
                print "One line, node1:(%f,%f), node2:(%f,%f), is too short ( < 1.0e-12 )." % (line.node1[0], line.node1[1],line.node2[0], line.node2[1])
                exit()
            self.minLineLen = min(self.minLineLen,line.length)

    def calc(self, threadNum=1):
        print " Start BEM calculation ".center(50, "=")
        # discrete 'self.freqs' List 
        numOfFreq = int(len(self.freqs)/threadNum)
        for ind in range(threadNum-1):
            exec("divedFreq%d = self.freqs[%d*numOfFreq:(%d+1)*numOfFreq]" % (ind, ind, ind))
        exec("divedFreq%d = self.freqs[%d*numOfFreq:]" % (threadNum-1, threadNum-1))

        # prepare threads.
        for ind in range(threadNum):
            exec("t%d = BemThread(self, divedFreq%d)" % (ind, ind))
            exec("t%d.start()" % ind)
        for ind in range(threadNum):
            exec("t%d.join()" % ind) # waiting for all thread finished.
        # Finish successfully
        print " Finish BEM calculation ".center(50, "=")

    


class BemThread(threading.Thread):
    """
    """
    def __init__(self, BEM, freqs):
        threading.Thread.__init__(self)
        self.freqs = freqs
        self.BEM = BEM

    def run(self):
        # freq loop start 
        for freq in self.freqs:   
            waveNum = 2.0*pi*freq/self.BEM.c # waveNumber
            waveLen = self.BEM.c/freq # waveLength
                
            # discrete boundary elements
            #-- maximum length of boundary elements, Ref [1] --
            limElmLen = min(waveLen/self.BEM.ndiv,2.0*self.BEM.minLineLen)
            #limElmLen = .01
            
            boundElm = self.discBoundElm(self.BEM.lines, limElmLen)
            print "%s : freq = %6.f[Hz], %d elements..." % (self.getName(), freq, len(boundElm))
            #BEM.alert("%s : freq = %6.f[Hz], %d elements..." % (self.getName(), freq, len(boundElm)))
            #-- error trap --
            if self.BEM.limElmNum < len( boundElm ):
                print "%s : Too many elements, no computation has done." % getName()
                #BEM.alert("%s : Too many elements, no computation has done." % getName())
                exit()
            
            # Build up and solve BEM calculation,and obtain velocity potentials on the elements. 
            solution = self.sbemEqn(self.BEM.srcPos, boundElm, waveNum)            
        
            for pos in self.BEM.recPos:
                posInd = find(((self.BEM.recPos[:,0] == pos[0]) * (self.BEM.recPos[:,1] == pos[1])) == True) 
                freqInd = find(self.BEM.freqs == freq)
                # Calculate SPL at the receiving points
                self.BEM.result[freqInd, posInd]=self.calField(self.BEM.srcPos, array(pos), boundElm, solution, waveNum)


    def discBoundElm(self, lines, limElmLen):
        """
        ======================================================
                 
         input>>
         Line : List of lineClass Instance (ref:boundClass.py)
         limElmLen : Maxmum limit of element length (>=1.0e-12 setted@sbemMain.py)
        
         output>>
         boundElm : List of elemntClass Instance
        
         Discrete boundaries to element
        
        ======================================================
        
         Function tree>>
         sbemMain -- discBoundElm
        
        ======================================================
        """
        # ================================ count Lines in form
        numLine=len(lines)
        # ================================ count the number pf elemnt & create memory space
        numElm=[0 for x in range(numLine)]
        for n in range(numLine):
            lineLen=linalg.norm(lines[n].node2-lines[n].node1) # scipyFunction : norm()
            numElm[n]=int(ceil(floor((float(lineLen)/float(limElmLen))*10**6)/10**6))
        # ================================ Discretise Lines
        # allocate memory
        boundElm=zeros((1,sum(numElm)))[0]
        boundElm=[element(0,0,0,0,0) for x in boundElm]
        # calc each position of nodes
        index=0 # total number of elements
        for n in range(numLine):
            k=linspace(0,1,numElm[n]+1) # From 0 to 1 evenly spaced numbers by numElm.
            tmpNode = vstack((lines[n].node1[0]+lines[n].vec[0]*k,
                              lines[n].node1[1]+lines[n].vec[1]*k))
            for m in range(numElm[n]):
                boundElm[index+m].node1=tmpNode.T[m]
                boundElm[index+m].node2=tmpNode.T[m+1]
                boundElm[index+m].admit=lines[n].admit
            index += numElm[n]
        return boundElm
    
    
    def sbemEqn(self, srcPos, boundElm, waveNum):
        """
        ======================================================
         def sbemEqn(srcPos, boundElm, waveNum)
        
         input>>
         srcPos : array_ positions of source sound. ref.)sbemMain.py
         boundElm : List_ classInstances of boundary Element. See also, boundClass.py & boundElm.py
         waveNum : global num_ A Wave Number.  ref.)sbemMain.py
         
         output>>
         solution : array_ velocity potential on the elements 
         
         Build up and solve BEM simultaneous equation.
         and velocity potentials at the middle point of boundary elementsare obtained.
        
        =======================================================
        
         Function tree>>
         BemThread --- sbemEqn---element.calcPara()
                         |
                         +-elmIntGen
                         +-elmIntRonElm
                         |
                         +-greenFunc
        
        =======================================================
        """    
        numElm=len(boundElm)
        # coefficient Matrix================================
        # prepare matrix
        # Dont't use "zeos". Because numpy's zeros can't receive complexNumbers.
        directVec=[[0.] for x in range(numElm)]#zeros((numElm,1))
        coeffMat=[[0. for x in range(numElm)] for x in range(numElm)]

        
        for col in range(numElm): # relation to ElementNumber
            boundElm[col].calcPara() # calculate each Parameter
    # loop for Element
        for row in range(numElm): # relation to ElementNumber        
            for col in range(numElm): # relation to ReceiverNumber
               if col != row:
                   coeffMat[row][col]=self.elmIntGen(boundElm[col],boundElm[row].mid,waveNum)               
               else: # col is row
                   coeffMat[row][col]=0.5+self.elmIntRonElm(boundElm[col],waveNum)               
            #greenFunc
            directVec[row]=self.greenFunc(srcPos,boundElm[row].mid,waveNum)
        
        solution=linalg.solve(coeffMat, directVec)    
        return solution # Return N*1 Vec
    
    
    def calField(self, srcPos, recPos, boundElm, solution, waveNum):
        """
        ======================================================
         def calField( srcPos, recPos, boundElm, solution, wavenum ):
        
         input>>>
         srcPos : [x,y]_ A list of source position
         recPos : [x,y]_ A list of receiver position ( may define plural receivers )
         boundElm : A list of boundClassInstances 
         solution : N*1 vec_ velosity potentials at the middel points of the elements
         wavenum : wave-number
        
         output>>>
         result : calculated normalized SPL
        
         Calculate normalized SPL at each receiver positions,
         using velosity potentials on the elements.
        
        ======================================================
        
         Function tree:
         BemThread --- calFiled -+- elmIntGen 
                                 |
                                 +- greenFunc
                    
        ======================================================
        """ 
        #=========================== 1m point potential in free field for normalization
    #    freeField1m=-0.25j*special.hankel2(0,waveNum)
        freeField1m=-0.25j*self.H20(waveNum)
        # loop for each receiving points
        scattered=0
        for n in range(len(boundElm)):
            scattered = scattered + self.elmIntGen(boundElm[n],recPos,waveNum)*solution[n]
        potential=self.greenFunc(srcPos,recPos,waveNum)-scattered
        result=20*log10(abs(potential/freeField1m))        
        return result
    
    
    
    def greenFunc(self, srcPos, reccPos, waveNum):
        """
        ======================================================
         
         input>>>
         srcPos : [x,y]_ A list of source position
         reccPos : [x,y]_ A list of receiver position ( may define plural receivers )
         wavenum : wave-number
        
         output>>> 
         vaule : GreenFunction's value
        
         Calculate a contribution from a line source floating in the free field
         to a receiver point.
        
        ======================================================
        
         Function tree:
         BemThread -+- sbemEqn --- greenFunc 
                    |
                    +- calField --- greenFunc
        
        ======================================================
        """
        pr=srcPos-reccPos # p-r : vector from ReceiverPoint to SpurcePosition
        kpr=waveNum*linalg.norm(pr) # k|p-r|
    #    value=-0.25j*special.hankel2(0,kpr)
        value=-0.25j*self.H20(kpr)
        return value
    
    
    def elmIntGen(self, rowElm, recPos, waveNum):
        """
        ======================================================
                 
         input>>
         rowElm[row] : boundClassInstance_ An instance in the boundElm.
         recPos : recording Position.
         waveNum: global num_ A Wave Number.  ref.)sbemMain.py
        
         output>>
         value : Local Num_ An integral value  
        
         calculate element integran value
         when the rec point is not located on the same element.
        
        ======================================================
        
         Function tree>>
         BemThread --- elmIntGen --- GLquad --- integrand
             |                   
             |
         calField --- elmIntGen --- GLquad --- integrant
                   
         
        ======================================================
        """

        #value = 0
        value = self.GLquad(-1,1, self.integrand, rowElm, recPos, waveNum) 
        value = 0.125*waveNum*rowElm.length*value
    
        return value
    
    
    def elmIntRonElm(self, colElm,waveNum):
        """
        ======================================================
                 
         input>>
         colElm[col] : boundClassInstance_ An instance in the boundElm.
         waveNum: global num_ A Wave Number.  ref.)sbemMain.py
        
         output>>
         value : Local Num_ An integral value  
        
         calculate element integrand value
         when the colElm.mid point is located on the same element.
        
        ======================================================
        
         Function tree>>
         sbemMain --- elmIntRonElm -+- integrant
             |
         calField -+- integrant
         
        ======================================================
        """
        # the element is not rigid
        if colElm.admit!=0:
            value=(colElm.admit*0.5)*self.singularIntH20(waveNum*colElm.length*0.5)
            # elmPara.admit == 0 , i.e. the element is rigid
        else:
            value=0
        return value
    
    def singularIntH20(self, a):
        """
        ========================================================== 
         value = singularIntH02( a )
         value : singular integral value 
         a : range parameter
         
         Function name : 'singularIntH02'
        
         Calculate singular integral value , \Int H02(z)dz ( 0 <= z <= a ).
        
        ========================================================== 
        
         Function tree:
         elmIntRonElm --- singularIntH02 
        
        ========================================================== 
        """
        var = 2j/pi # 2j/pi
        ha = 0.5*a # a/2
        logha = log(ha) # log(a/2)
        #======================================== Sum terms ( m <= 20 )    
        # m=0 =================
        xim=-0.57721566 # - Euler's constant 
        factm=1 # m!
        sign=1 # (-1)^m
        value= 1 - var*(logha-xim-1)
        # 0<m<=20 =============
        for m in range(1,21):
           xim = xim + (1.0/m)
           factm = factm * m
           sign = sign * (-1)       
           termF = ( sign*(ha**(2*m)) ) / ( (factm**2) * (2*m+1) )
           termS = logha - xim - ( 1/(2*m+1) )
           value = value + ( termF * ( 1 - var * termS ) )
        
        value=value*a;
        return value
    
    def GLquad(self, low, up, func, *var):
        """
        ======================================================
                
         input>>
         low : lower limit of integration.
         up : upper limit of integration.
         func : An object refers to function.
         *var : Variable numbers for function.
        
         output>>
         value : A result of integration of func from "low" to "up" 
        
        
        ======================================================
        
         Function tree>>
         elmIntGen --- GLquad --- integrand
         
        ====================================================== 
        """ 
         
        u=[-0.932469514203152,-0.661209386466265, -0.238619186083197, 
            0.238619186083197, 0.661209386466265, 0.932469514203152]
        w=[0.1711324492379170, 0.360761573048139,0.467913934572691,
            0.467913934572691, 0.360761573048139, 0.1711324492379170]
        
        #Calculate GL integration 
        value=0
        divWidth = (up - low)/2.0
        midWIdth= (up+low)/2.0
        for  i in range(len(u)):
            x = divWidth*u[i] + midWIdth
            value += w[i]*func(x, *var)
        value= value*divWidth
        return value
    
    def integrand(self, t,rowElm,recPos,waveNum):
        """
        ======================================================
                 
         input>>
         mode : change the function by the boundElm.admit. For reduce calculation volume.
         t : num_ A number of integer variable Number.
         rowElm : A boundClass Instance.
         recPos : [x,y]_ A list of receiver Point.
         waveNum : 
        
         output>>
         boundFunc : functionReference_  Function suffering integral calculus 
        
        
        ======================================================
        
         Function tree>>
         elmIntGen --- GLquad --- integrand
         
        ====================================================== 
        """   
        # prepare some numbers. Please check out "BEM method"
        qr=rowElm.mid+rowElm.range*t - recPos # q-r
        absQr=linalg.norm(qr) # |q-r|
        kqr=waveNum*absQr # k|q-r|
        ipr=dot(rowElm.normalVec,qr) # innerProduct_ n (dot) q-r
        dcos=ipr/absQr # Directional cosine
        # make Function
        if rowElm.admit is not 0: # bound is not rigid. rowElm.admit!=0 
    #        intFunc=1j*dcos*special.hankel2(1,kqr)+rowElm.admit*special.hankel2(0,kqr)
            intFunc=1j*dcos*self.H21(kqr)+rowElm.admit*self.H20(kqr)
        else: # bound is rigid. rowElm.admit=0
    #        intFunc=1j*dcos*special.hankel2(1,kqr)
            intFunc=1j*dcos*self.H21(kqr)
        return intFunc
    
    
    
    def H10(self, x):
        """
        ======================================================
        HankelFunction of the first kind of oder 0
        For the use of reducing time to calculate.
        H10(x) ~= hankel1(0, x)
        
        0<x<3.0 : Ref[1] p. 369 9.4.1, 9.4.2 \n
        x>=3.0  : Ref[1] p. 369 9.4.3 \n
        
        Computation Error:
        |x|<3.0 -> |Re{err}|<5.0e-8, |Im{err}|<1.4e-8
        |x|>=3.0 -> |Re{err}|<1.6e-8, |Im{err}|<1.6e-8
        
        If argument x is negative value, return 0.
        Other input, return "value"
        
        Reference:
        [1]:Milton Abramowitz and Irene A.Stegun, 
        "Handbook of Mathematical functions," Dover Publications, 1965.
        ======================================================
        """
        value = 0
       
        if x<=0: # x<=0         
            pass
        elif  3<=x : # 3<=x 
            px = x
            s = 3/px
            f = ((((( .00014476*s -.00072805 )*s +.00137237 )*s -.00009512 )*s-.00552740 )*s -.00000077 )*s +.79788456
            t = ((((( .00013558*s -.00029333 )*s -.00054125 )*s +.00262573 )*s -.00003954 )*s -.04166397 )*s -.78539816 + px
            value = ( math.cos(t) + 1j*math.sin(t) ) * f/math.sqrt(px) 
        else: # 0<x<3  
            px = x
            s = ( px * .33333333 ) 
            s = s*s
            j0 = ((((( .00021*s -.0039444 )*s +.0444479 )*s -.3163866 )*s +1.2656208 )*s -2.2499997 )*s + 1.
            y0 = ((((( -.00024846*s +.00427916 )*s -.04261214 )*s +.25300117 )*s -.74350384 )*s +.60559366 )*s +.36746691 + j0 * math.log(px*.5)*(2/math.pi)
            value = j0 + 1j*y0
        return value
    
    
    def H20(self, x):
        """
        ======================================================
        HankelFunction of the second kind of oder 0
        For the use of reducing time to calculate.
        H20(x) ~= hankel2(0, x)
        
        0<x<3.0 : Ref[1] p. 369 9.4.1, 9.4.2 \n
        x>=3.0  : Ref[1] p. 369 9.4.3 \n
        
        Computation Error:
        |x|<3.0 -> |Re{err}|<5.0e-8, |Im{err}|<1.4e-8
        |x|>=3.0 -> |Re{err}|<1.6e-8, |Im{err}|<1.6e-8
        
        If argument x is negative value, return 0.
        Other input, return "value"
        
        Reference:
        [1]:Milton Abramowitz and Irene A.Stegun, 
        "Handbook of Mathematical functions," Dover Publications, 1965.
        ======================================================
        """
        
        value = 0
       
        if x<=0: # x<=0         
            pass
        elif  3<=x : # 3<=x 
            px = x
            s = 3/px
            f = ((((( .00014476*s -.00072805 )*s +.00137237 )*s -.00009512 )*s-.00552740 )*s -.00000077 )*s +.79788456
            t = ((((( .00013558*s -.00029333 )*s -.00054125 )*s +.00262573 )*s -.00003954 )*s -.04166397 )*s -.78539816 + px
            value = ( math.cos(t) - 1j*math.sin(t) ) * f/math.sqrt(px)
        
        else: # 0<x<3  
            px = x
            s = ( px * .33333333 ) 
            s = s*s
            j0 = ((((( .00021*s -.0039444 )*s +.0444479 )*s -.3163866 )*s +1.2656208 )*s -2.2499997 )*s + 1.
            y0 = ((((( -.00024846*s +.00427916 )*s -.04261214 )*s +.25300117 )*s -.74350384 )*s +.60559366 )*s +.36746691 + j0 * math.log(px*.5)*(2/math.pi)
            value = j0 - 1j*y0
        return value
    
    #===================================================
    
    def H21(self, x):
        """
        ======================================================
        Hankel function of the second kind of oder 1
        H21(x) ~= hankel2(1,x)
        
        0<x<3.0 : Ref[1] p. 369 9.4.1, 9.4.2 \n
        x>=3.0  : Ref[1] p. 369 9.4.3 \n
        
        Computation Error:
        |x|<3.0 -> |Re{err}|<1.3e-8, |Im{err}|<1.1e-7
        |x|>=3.0 -> |Re{err}|<4.0e-8, |Im{err}|<4.0e-8
        
        If argument x is negative value, return 0.
        
        Reference:
        [1]:Milton Abramowitz and Irene A.Stegun, 
        "Handbook of Mathematical functions," Dover Publications, 1965.
        ======================================================
        """
        value = 0
       
        if x<=0: # x<=0         
            pass
        elif 3<=x: # 3<=x 
            px = x
            s = 3 / px
            f = ((((( -.00020033*s +.00113653 )*s -.00249511 )*s +.00017105)*s +.01659667 )*s +.00000156 )*s +.79788456
            t = ((((( -.00029166*s +.00079824 )*s +.00074348 )*s -.00637879)*s +.0000565 )*s +.12499612)*s -2.35619449 + px
            value = ( math.cos(t) - 1j * math.sin(t) ) * f / math.sqrt(px)        
        else: # 0<x<3
            px = x
            s = ( px * .333333333 )
            s = s * s
            j1 = (((((( .00001109*s -.00031761 )*s +.00443319 )*s -.03954289)*s +.21093573 )*s -.56249985 )*s +.5 ) * px
            y1 = (((((( .0027873*s -.0400976 )*s +.3123951 )*s -1.3164827 )*s +2.1682709 )*s +.2212091 )*s -.6366198 ) / px + j1 * math.log(px*.5) * (2/math.pi)
            value = j1 - 1j*y1
        return value
    
def attributesFromDict(dict):
    """set attributes to this object.
    this code refer to 'PythonCookBook' p284 from O'reilly in Japan.
    """
    self = dict.pop("self")
    for name,value in dict.iteritems():
        setattr(self, name,value)    
# [e.o.f]